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Abstract 
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pg I The generation of intrinsic rotation by turbulence and neoclassical effects in tokamaks is con- 

sidered. To obtain the complex dependences observed in experiments, it is necessary to have a 
^7^1 model of the radial flux of momentum that redistributes the momentum within the tokamak in 



c/3 
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the absence of a preexisting velocity. When the lowest order gyrokinetic formulation is used, a 
symmetry of the model precludes this possibility, making small effects in the gyroradius over scale 



c/2 . length expansion necessary. These effects that are usually small become important for momentum 

c/3 , transport because the symmetry of the lowest order gyrokinetic formulation leads to the cancella- 

"^ ' tion of the lowest order momentum flux. The accuracy to which the gyrokinetic equation needs to 



be obtained to retain all the physically relevant effects is discussed. 
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I. INTRODUCTION 



Tokamak plasmas can rotate toroidally to high speeds because they are axisymmetric. 
The presence of toroidal rotation is beneficial: it suppresses large-scale MHD instabilities 
pj, |2| and reduces the turbulent radial fiux of particles^ energy and momentum 3|-l8J. In 
a number of experiments (e.g. JET J9| and Dlll-D [10]) toroidal rotation is achieved by 
direct injection of momentum via neutral beams. In these cases the rotation can be close 
to sonic, resulting in turbulence suppression by velocity shear and a transition to regimes of 



reduced transport 



m 



12| . However, neutral beam injection of momentum is less effective 



in large, dense tokamak plasmas, and in particular, it is expected that the momentum 
input in ITER [13] will be small. This situation has driven a great deal of research on 
spontaneous, or intrinsic, rotation [14]. Even in experiments with no momentum injection, 
it is often observed that the toroidal rotation is non-negligible. This intrinsic rotation shows 
non-trivial dependences on different plasma parameters. 

Intrinsic rotation can be in the direction of the current or against it, and often both signs 
are found in the same shots at the same time in different locations within the plasma 15l. Il6[. 
The rotation at the edge is usually co-current and it is not correlated to the rotation in the 
core, i.e. it is possible to find pulses with almost identical rotation at the edge that will 
have opposite velocities in the core. Several experimental measurements indicate that the 
gradients in temperature and pressure are the drive for intrinsic rotation 17H19I]. However, 



it is found that pulses with very similar density and temperature profiles have intrinsic 
rotation profiles that have the opposite sign, and it is possible to jump from co-current to 
counter-current rotation within the same pulse with slight changes in the density and the 
temperature. These rotation reversals seem to be related to a density threshold 20|, |2lj in 



TCV 



22| and Alcator C-Mod [23] , but it is not clear what determines the direction of the 



rotation in other tokamaks such as JET 2J|. RF heating and current drive complicate the 
picture even more because Ion Cyclotron Resonance Heating (ICRH) 25l-l27l] tends to drive 
co-current rotati on, and Electron Cyclotron Resonance Heating (ECRH) drives counter- 
current rotation llSl^Sl. Lower Hybrid (LH) current drive can drive co- or counter-current 

rinn 

rotation depending on the situation [26|, |29|, |30|]. Any modeling effort needs to consider 
several effects and must allow them to compete in order to obtain the variety of intrinsic 
rotation profiles observed in experiments. 



In this article we discuss the features required for a self-consistent model of intrinsic ro- 
tation. Since rotation is determined by the radial flux of toroidal angular momentum, any 
model of intrinsic rotation will necessarily have to include the main cause of radial trans- 
port: the turbulence. The need to model the turbulence makes the problem particularly 
challenging because turbulence is a rich and complex phenomenon, and in tokamak plasmas 
is kinetic in nature. The problem is simplified by employing gyrokinetics [3l|, a reduced 
kinetic model adequate for low frequency fluctuations. Gyrokinetics averages over the gy- 
rofrequency, which is too fast to matter for the fluctuations, and at the same time keeps 
wavelengths that are comparable to the gyroradius. To do this, an expansion in the small 
parameter p^, = pi/L '^ 1 is performed, and terms of higher order neglected. Here pi is the 
ion gyroradius and L is the typical size of the machine of interest. 

The problem with modeling intrinsic rotation is that the momentum flux that determines 
the intrinsic rotation profiles can only be obtained by solving for small corrections to the 
lowest order gyrokinetic treatment. This is necessary because of a symmetry present in 
the lowest order gyrokinetic system, first proposed in the context of quasilinear estimates 



for turbulent transport 32[. The complete symmetry of the non-linear set of gyrokinetic 



equations has only been recently derived and confirmed with gyrokinetic simulations 331.134 1. 
Because this symmetry was only recently identified, the turbulence community had not 
worried about the order to which the expansion in p^ had to be carried out to obtain the 
correct radial flux of toroidal angular momentum. The order in p* to which is necessary to 



solve the gyrokinetic equation was obtained in [35||, where it was assumed that there is a 
cancellation in the radial flux of momentum due to up-down symmetry. 

One of the main conclusions of 35|] is that none of the existing simulations of turbu- 
lent momentum transport contain all the physical phenomena that can generate intrinsic 
rotation. How far the codes are from a comprehensive and complete model depends on 
the approach. Some codes use fluid equations to evolve the proflles of density, temperature 



and rotation 
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371], and the gyrokinetic equation is solved in a 6f code |38l-l42| only to 
obtain the self-consistent fluctuations to then calculate the radial fluxes of particles, energy 
and toroidal angular momentum. In these simulations both the fluid equation for toroidal 
angular momentum and the gyrokinetic equation solved in the 6f code have to be modifled 
to include effects that are an order higher in p*. The problem is different for simulations 
that solve for the full distribution function and do not use fluid equations to evolve the 



density, the temperature and the rotation |43|-|46| . In these codes the rotation profile is 



obtained from a moment of the distribution function, and the transport equations that were 
used in the other approach are assumed to be contained correctly by the gyrokinetic Fokker- 
Planck equation for the full distribution function. It is far from obvious that the gyrokinetic 
equation describes correctly the evolution of density, temperature and velocity because the 
equation for the full distribution function is a gyrokinetic equation and it is truncated to 



some low order in p*. This limitation was pointed out in a series of articles 47H50| . and 
it was concluded that it is necessary to keep many more terms in the gyrokinetic equation 
than the terms retained in existing simulations to recover the correct conservation equation 
for the toroidal angular momentum. 

In this article we will review all these results and explain how they apply to the modeling 
of intrinsic rotation. In Section |TT] we will discuss the approach that uses fluid equations 
because it is the most intuitive path to obtain the order to which the equations have to be 
expanded. In Section |TTT] we discuss the requirements on the approaches that do not use 
fluid equations, and in particular, we derive the implications for collisionless models that 
employ a phase-space Lagrangian gyrokinetic formalism. Finally, in Section IIVI we close 
with a discussion that includes possible future approaches to the problem. 

II. INTRINSIC ROTATION 

In most tokamaks the time between ion-ion collisions is shorter than the transport time 



scale. As a result 
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48|, the ion velocity is given by the neoclassical formula 51 1 



where the toroidal rotation angular velocity is 

and fii is a numerical factor that depends on the geometry of the flux surface and the 
collisionality Rua/vu- Here n^, Vj, Tj and Pi = riiTi are the ion density, average velocity. 



temperature and pressure, 0o is the lowest order electrostatic potential, vu = j2Ti/mi is 
the ion thermal speed, va is the ion-ion collision frequency, Ze and mi are the ion charge 
and mass, R is the major radius, C, is the unit vector in the toroidal direction, %p is the 



poloidal magnetic flux and has the role of radial coordinate, ( is the toroidal angle, B = 
JVC + VC X Vip and B are the magnetic field and its magnitude, / = -R-B,^, -B,^ is the 
toroidal component of the magnetic field, VC = C/Ry ^ind e and c are the proton charge and 
the speed of light. To obtain the rotation fi^ we need to solve the conservation of toroidal 
angular momentum, given by 

|(y'(i?n.m.V..C)J=-A(I/'n) + 5,. (3) 

Here (. . .)^ = {Vy^ J dO dC {B ■ Ve)-\. . .) is the flux surface average, V = dV/diP = 
J d9 dC (B • V9)~^ is the volume contained between two contiguous flux surfaces, 9 is the 
poloidal angle, 5*^ is a momentum source and 



n = /(i?m,|d3t;/,(vO(vV^) 



-i-i?(B.C)(B-V^))^)^ (4) 

is the radial flux of toroidal angular momentum, where /j(r, v, t) is the ion distribution 
function, and (• • .)t = (^^)"^(^V^) ""^ /a* '^^/ai/. '^'^(- • •) i^ the coarse grain average over a 
time interval At and a radial segment Aip to average over the turbulent fluctuations. The 
radial flux of toroidal angular momentum has two main components: the motion of parti- 
cles that transfers toroidal angular momentum from one flux surface to the next, given by 
Rrrii j d?v fiiy • C,){y ■ VV"), and the Maxwell stress due to the fluctuating magnetic field 
(47r)~^i?(B ■ ^)(B • V^) (note that the non- fluctuating part of the magnetic field satisfies 
(B)t ■ V'?/' = 0). In general, for low /3 plasmas the Maxwell stress contribution is small 
compared to the radial flux of momentum due to the particles. For the rest of this article we 
will neglect it. We will also assume that the electric field is electrostatic E = — V0. These 
assumptions are easily relaxed and do not change the final result. 

Even in the absence of momentum injection, it is observed that tokamaks have a sizable 
rotation. To study this intrinsic rotation we need to find a non- zero solution to ([3]) for S*,^ = 0. 
If in addition we consider a steady state, we are left with the equation d{y'I\)/dil) = that 
implies that V'U is a constant. By imposing regularity at the magnetic axis, where Vip = 0, 
the final equation for the intrinsic rotation is that there is no net momentum flux through 
any flux surface, i.e. 

n[(]c(V^);ne(^),T,(^),T,(^),...]=0. (5) 



Thus, to solve for the profile of intrinsic rotation we need to obtain the functional dependence 
of n on ^(^{ip) and then solve the equation 11 = 0. 

The radial flux of momentum 11 is dominated by turbulence. In the following sections 
we will treat the turbulence by expanding in the small parameter p^, = Pi/L ^1. In 
subsection III Al we review the characteristics of the turbulent fluctuations in a tokamak and 
we use it to split the ion distribution function and the electrostatic potential into different 
pieces. In subsection III Bl we write the radial flux of toroidal angular momentum in a specific 
form that in combination with the ordering of subsection III Al will give a complete picture 
of the necessary pieces for momentum transport. We will find that one of the pieces in the 
radial flux of momentum seems to dominate. However, in subsection III CI we are able to 
show using a symmetry of the system that the apparently large term in the radial flux of 
momentum vanishes and we need to keep the other terms. We discuss the implications for 
intrinsic rotation in subsection IIIDI 

A. Ordering 

In tokamaks the turbulent fluctuations are a factor of p=K ^ 1 smaller than the background 
quantities. Their characteristic frequency is of the order of vu/L, and their characteristic 
wavelength is of the order of the ion gyroradius pj. The scale length of the background 
quantities is much longer, of the order of the size of the device L, and they evolve on a 
much longer transport time scale te ~ p~'^L/vti ^ L/vu- Thus, it is useful to think about 
the distribution function as composed of a long wavelength piece Fi and a short wavelength 
turbulent component ff", i.e. 

/, = F,(V',^,v,t) + /f(r,v,t). (6) 

Similarly, the electrostatic potential can be split into a long wavelength piece $ and a 
fluctuating piece 0*^, leading to 

= $(^,^,t) + 0*^(r,t). (7) 

The long wavelength piece of the distribution function that contains the background 
density, temperature and velocity is axisymmetric {d/dC, = 0) and is composed of pieces of 
different order in p*, giving 

F^ = fM^ + i^r + F.2 + ..., (8) 



where fMi is the lowest order MaxweUian, F^^ ~ p^fui is the next order correction that 
contains the neoclassical corrections |51i, and Fj2 ~ plfhii is the second order correction that 
has to do with both neoclassical processes and the turbulence. In |35| these different pieces 
are calculated employing a simplifying expansion in the ratio of the poloidal component of 
the magnetic field Bp and the total magnetic field 5, B^/ B <^ 1. The equations for the 
complete Fj2 without the expansion in Bp/ B <^ 1 are given for the first time in [521 ] . 

The turbulent pieces of the potential and the distribution function have perpendicular 
wavelengths comparable to the ion gyroradius, k±pi ~ 1, but the parallel gradients are much 
smaller, k\\L ~ 1. This anisotropy is due to the critical balance between the characteristic 
frequency of the parallel motion k\\Vti and the characteristic frequency of the perpendicular 
drifts fc_LV^ ~ kj_p^Vti, where v^ = —{c/B)'V(j)^^ x b ~ p*Vti is the turbulent E x B drift 
53| . The fact that the characteristic perpendicular size of the eddies is much smaller than 
the characteristic size of the tokamak suggests the use of a local approximation in which the 
fluctuations are treated as an eikonal, i.e. 

0*'^(r, t) = J2 0*''(^, 0, k^, fc„, t) exp(?k^ ■ r) (9) 

and 

/*"(r, V, t) = ^ f^'\i:, e, k^, k^, V, t) exp(2kx ■ r), (10) 

kx 

where the perpendicular wavenumber is 

k^ = k^Vip + kaVa . (11) 

pi 

Here i = a/— 1, and a is a magnetic coordinate perpendicular to the magnetic field and 

parallel to the flux surface defined by B = Va x Vip. The fluctuations have a fast spatial 

variation associated with k±_ and a much slower poloidal and radial variation. The slow 

poloidal variation is obtained for ip and a fixed, and has to do with the long correlation 

lengths along magnetic field lines. The fluctuations (p^^ and f^^ can be expanded in p* as 

well, giving 

0*' = 0f + <pf + ... (12) 

and 

/^' = /A' + /^V..., (13) 

where ec^f/T, ~ /^f //m. ~ p. and e^^V^e ~ H/fNU ~ pl 



B. Radial flux of toroidal angular momentum 

We will see soon that the radial flux of toroidal angular momentum 11 has to be calculated 



up to order plpiR\W'4)\ to describe intrinsic rotation. In 35| a convenient form of the radial 



flux of toroidal angular momentum is calculated using moments of the Fokker-Planck equa- 



tion. The convenience of this expression is discussed in Appendix B of [5J]. To summarize 
the result found in [35,] , the radial flux of momentum is made of two pieces 

n = n_i + Ho. (14) 

The piece n_i is formally of order p1piR\Wip\ and is given by 

n-i = ^^m,ci?^/d\/AV-C)^ 



^ " I Ip I 'Y 



J^'vC\^[F:^]{v.C,?) , (15) 



2Ze . , , 



where C^-j [/] is the linearized ion-ion collision operator. The contribution Ho is of order 

P*Pi-R|V'?/'| and is given by 



Ii, = lUcR^-^jd\m^-C) 



dC 

Lth 



ipl X 



Ipl T 



-(^/d3.Cr[..](v.C)^)^ 



' mfcR^ 



j^vC^'^if^jn^-CY 



\\ 2Ze . I ,, 

\ \ / ill r^ 

1 d ^^J lm'ic^R^d<pf " ^ ^ 



+T^5^^'Ulz^^y'^'"^"(-«'^ 



where C|" [/,(?] is the complete bilinear ion-ion collision operator. If n_i does not vanish, 
it dominates and it is unnecessary to compute the piece IIq and hence the second order 
corrections Fj2, f^ and 02*^. However, we will show in the next subsection that n_i vanishes 
for most tokamaks when momentum is not injected into the plasma due to a lowest order 
symmetry of the turbulence and the neoclassical transport. 

8 



C. Up-down symmetry of the flux of toroidal angular momentum 



The lowest order neoclassical and gyrokinetic equations determine the functions 
Fl[^(tp,9,v\\, fi,t), f^._^{'ilj,9,k^,ka,v\\, ^,t) and 0* {il),9,k^,ka-,t), where the velocity space 
is described by the parallel velocity v\\ and the magnetic moment /i = v\/2B. These func- 
tions depend on the local density and temperature gradients, the local rotation f2^, and its 
gradient dVl^/dil) at the location ^. The functions Ff^{il}, 6, v^, fi,t), /^^{"ip, d, k^, ka, v\\, ^, t) 
and 0* {ip,9,k^,ka,t) only depend on the local parameters because the characteristic size 
of the turbulence and the width of the drift orbits that lead to neoclassical transport are 
small compared to the size of the machine. 

The rotation and the rotation gradient appear in the lowest order equations for 



F/}=(?/^,6',W||,/i,t), l.^{^,e,k^,ka,vii,iJ,,t) and 0^ {ilj,e,k^,ka,t) onf 
dered to be sonic, RQr ~ Vu- This limit is the high flow ordering [55 



the rotation is or- 
, and it has been 



adopted by many of the numerical studies of turbulent momentum transport performed so 



far 



59 



62| . If the rotation is sonic, most of the terms in ([T]) and (|2]) can be neglected 
to lowest order in p*, leaving only Vj = Q(^R(^, with Q(^ = —c{d(j)o/d4'). When the high flow 
ordering is employed, the lowest order equations for F^{il), 6, f ||, /i, t), f\^ {ip, 9, k^, k^, f ||, /i, t) 
and 0* {ip, 9, k^, ka, t) are invariant under the following transformation in an up-down sym- 
metric flux surface: the sign of the rotation fi^ and the rotation shear dQ(^/dtp is reversed, 
the sign of the variables f ||, /i and k^ is reversed, and the sign of the functions F^l^, /*^ and 



,tb 



is reversed, i.e. 



Fl{'{ij,9,v\\,fi,t;n^,dnjd^)-> 

-Fii^ij, -9, -v\\,fi, t- -n^, -dnjdij), 



(17) 



and 



/.^ (^, 9, /c^, A;„, vii, /x, t; n^;, dnjd^) -^ 



-Lii (^' ~^' -h^ ^"' -^11, yw, t; -fif , -dnjd^/j) 



^^{ilj,9,k^,ka,t-^od^dd^) 



-)■ 



0f (t/;, -9, -k^, ka, t; -n^, -dnjd^). 



discussed in 



(18) 



(19) 



Symmetry flT7|) can be found in Appendix A of [35], and symmetries (fT8l) and (TT9!) were 



34| . When these symmetries are applied to the lowest order piece of the 



moinentum transport n_i in (fTSj) . we find that in up-down symmetric flux surfaces, n_i is 
odd in Q^ and dVLc^/dip, 

n_i(fi^, dvijdti)) = -n_i(-fi^, -dnjdip), (20) 

leading to 

n_i(f]^ = 0, onjd^ = 0) = 0. (21) 



Tliis symmetry is broken by up-down asymmetry in tlie flux surface 63|, |6J], a fact that 



has been checked experimentally 65|]. The effect of asymmetry is small in most existing 
tokamaks, particularly in the core, where the magnetic flux surfaces are very up-down sym- 
metric. 

D. Consequences for intrinsic rotation 

Equation (12T|) shows that there is no momentum flux in the absence of rotation, that is, 
without preexisting rotation, the turbulence or neoclassical effects will not redistribute the 
momentum. Another way of looking at this is that a solution for n_i(i7(^, dQ(^/dip) = with 
the boundary condition i7^ = at the edge is simply Q(^{ip) = everywhere. There is no 
preferred direction and the plasma does not know in which direction it should rotate. This 
lack of preferred direction is only true for the high flow ordering. When we allow for a more 
general ordering below, we will be able to redistribute momentum even when there is no 
initial rotation. 

If the boundary condition is not ^2^ = at the edge, the solution to n_i(r2(^, dQt^/dtp) = 
can be some non-trivial rotation. However, the proflles allowed by a flux of momentum odd 
in rotation and rotation shear are very limited, and will only depend on the condition at the 
edge. To see this, we Taylor expand around fi,^ = and dVL^/dip = to obtain 

n-i ^ -P^c -x^ = o. (22) 

Note that physically the high flow ordering only has diffusion and a pinch J59l] in up-down 
symmetric configurations. Solutions to this equation are always of the form 

ficWocexp(^-/d^^). (23) 

This rotation profile never changes sign from co-current to counter-current, making the high 
flow model unable to predict many observed proflles that change sign [l5|, [l6|. Another 
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undesirable feature of solutions of the type (123|) is that they depend completely on the 
boundary condition at the edge. This contradicts observations such as the rotation reversals 



20 



2l| in which a small change in the density leads to a jump from a completely co-current 
profile to a profile in which the rotation at the magnetic axis is counter-current and the 
rotation at the edge remains the same, that is, co-current. 

To summarize, the high fiow ordering, for which i?f2^ ~ Vu, gives n_i 7^ and hence 
makes this piece of momentum transport the dominant contribution. However, it is not pos- 
sible to explain the complex dependences observed in experiments with the results obtained 
with the high fiow ordering. The reason is that the high flow ordering in up-down symmetric 
configurations only contains diffusion and a pinch, and therefore lacks any redistribution of 
momentum that is not induced by the presence of preexisting velocity. 

The inability of the high flow ordering to cope with the intrinsic rotation problem indicates 
that we need to allow the velocity to be RQ(^ ~ p*Vti. With this low flow ordering, O^ 
and dfi(^/dip do not enter in the equations for the lowest order pieces F^'^{ip,9,v\\, fi,t), 



f.^{ijj^6^k^,ka,v\\,^,t) and (j)^ {ijj,6,k^,k^t), and n_i is identically zero. As a result. Ho 



must be calculated. The formalism of 



33| proves elegantly that Ho is unaffected by the 



symmetry that made n_i = 0, and a contribution to the momentum transport. Hint, is 



independent of rotation 5J], giving 



afif , „ (Ale m m \ „,, 

Note that the solution to Ho = clearly gives non-trivial solutions for f2^ (■?/;). 

Thus, to correctly model intrinsic rotation, we need to obtain the contribution Ho ~ 
pl^PiR\Vip\, and this implies obtaining the second order pieces Fj2, //I" and 02*^. Several 
authors have made some attempts at calculating this higher order contribution to momentum 



transport. In 35|,|5J] a systematic expansion based on a subsidiary expansion in BjJ B ^ 1 is 
carried out. This subsidiary expansion assumes that the turbulent eddies and the size of the 
fluctuations does not scale strongly with B/Bp, and this simplifles the higher order equations 
enormously. Other authors have assumed that the main effect that contributes to Hnis the 



modiflcation of the turbulence by the radial variation of the background gradients 66N68|. 
Keeping only the radial variation of the background gradients means neglecting many other 
different effects that are as important as the one retained. Only by treating all the relevant 
terms on equal footing will it be possible to correctly reproduce the experimentally observed 
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dependences. An extension of the work done in 35|, |5J] that relaxes the assumptions on the 



turbulence but still relies on the Bp/ B ^ 1 expansion will be published shortly. This new 



35 



and in addition will 



improved model will include the effects already discovered in 

have the radial variation of the gradients and other effects that have not been considered 

before. 

III. MOMENTUM TRANSPORT IN GYROKINETICS 

From the results in the previous section, it is clear that to correctly predict intrinsic 
rotation, the radial flux of toroidal angular momentum must be calculated up to order 
p2pii?|V'?/'|. To do so, it is necessary to evaluate the next order correction to the fluctua- 
tions f^2 and (j)^2- These corrections are small compared to what is usually calculated in 
gyrokinetic simulations by one factor of p* ^ 1. 

A transport simulation that is based on fluid equations to evolve profiles can use the radial 
flux of toroidal angular momentum in f ll6p and then evaluate the turbulent fluctuations up 
to second order. To obtain the second order pieces, current gyrokinetic simulations have 



to be modified to include higher order effects such as those calculated in 35|, |5J]. These 
improvements are difficult to implement, but still tractable. 

Transport codes that evolve the full distribution function and then obtain density, tem- 
perature and velocity profiles from moments of the distributions function require much more 
work to recover the correct intrinsic rotation profile. In [49J the order of magnitude esti- 
mate n ~ p1piR\Vip\ was used to argue that these types of simulations need to obtain the 
gyrokinetic equation to order p'^fiVu/L, that is, to three orders of magnitude higher than it 
is usually calculated (only terms of order p^fiVu/L are retained in most gyrokinetic full / 



49|), but we can give an approx- 



simulations). The reasons for this estimate are subtle (see 

imate idea here. The electrostatic potential is calculated using a gyrokinetic quasineutrality 

equation that naturally contains a polarization density, 

V ■ (^ Vx0) +... = -Z jd?vUg + j d\ /e,. (25) 

In this equation we have kept a schematic polarization density on the left side, and the 
integrals over the distribution functions of the guiding centers of electrons, /eg, and ions, 
fig, on the right side. The distribution functions of guiding centers are evolved using gy- 
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rokinetic equations. This equation is only meant to be schematic and it is not intended to 
be exhaustive. It will only serve for an order of magnitude estimate. In full / simulations, 
the long wavelength electrostatic potential 0o is to be obtained using equation 1^5^. From 
([T]) and 02]) we find that Vj ~ {c/B)'V±(j)o. Then, the time derivative of the polarization 
density due to the background potential 0o is given by 



d_ 
di 



' f Zcrii „ 
V- ^Vj 

V5a 



L\eB'^' dt 



(26) 



and using the momentum transport equation ([3]) to make the order of magnitude estimate 
nimi{dVi/dt) ~ R~^dl[/dip ~ plpi/L, we finally find 



d_ 

di 



„ / Zcrii „ 



\Ba 

Using this estimate in (l25l) leads to 



''"'-t 



C 3Pi\ 4 ^ti 



i3 ^/eg 4 '^ti 



(27) 



(28) 



This condition implies that the gyrokinetic equations that give dfig/dt and dfeg/dt must be 
obtained up to order p^fiVu/L. This requirement makes the use of codes that are not based 
on fiuid equations to model intrinsic rotation an almost impossible task. 

The method that is most widely used to derive the gyrokinetic equation for full / sim- 
ulations is the phase-space Lagrangian approach [69|. For this reason, we now proceed to 
discuss the implications of what we have derived so far for Lagrangian methods. In the 
Lagrangian approach to the derivation of gyrokinetics, a coordinate transformation is per- 
formed in the phase-space Lagrangian of a particle in an electromagnetic field. The idea is 
obtaining a phase-space Lagrangian independent of gyrophase. This cannot be done exactly, 
but the gyrophase dependence can be eliminated to any order of interest (although the work 
required increases rapidly with the orders). Once the phase-space Lagrangian is calculated, 
a variational approach 70|, lDJ is used to obtain the equations of motion for the particles 
and Maxwell's equations for the electromagnetic fields. This variational method makes it 
necessary to have a Hamiltonian correct to order plT^ to obtain a quasineutrality equation 
sufficiently accurate to obtain /*^ and (p^^. The need for a Hamiltonian of this order to 



calculate the second order pieces fl2 and (p^^ has been argued for a slab in [49|] and for a 
tokamak in 1721. 



The question of the order to which the Hamiltonian is needed is not completely resolved 
by this argument. If the gyrokinetic equation is used only to calculate the fiuctuations. 
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which are then employed to evaluate the radial turbulent flux of momentum (fT6|) , we would 
not need to consider anything further. However, if the variational formulation is used to 
evolve the full distribution function and then the velocity is evaluated by taking a moment, 
we need to check if the gyrokinetic equation obtained with the Hamiltonian up to order plT^ 
is sufficient to recover the correct transport of momentum. The answer is not obvious since 
the fourth order Hamiltonian gives corrections to the gyrokinetic equation that are of order 



p^fiVti/L, and these may matter for momentum transport, as we saw in fl28|) . In 49|] it was 
proven that the Hamiltonian up to order plTe was sufficient in a slab. It is possible to show 
that it is also sufficient for a tokamak. By using the variational Lagrangian formulations, it 
is possible to obtain conservation equations that look similar to the conservation of toroidal 
angular momentum [73|, |7^ . The actual physical conservation of toroidal angular momentum 
and these gyrokinetic conservation equations are not equal to each other, but they tend to 
the same equation at long wavelengths. This property is what makes it possible to prove that 
the Hamiltonian to order plT^ is sufficient to obtain the correct rotation. The calculation 
that has achieved the highest order Hamiltonian is [72] , and the Hamiltonian in this reference 
is only accurate to order plT^. Thus, it is still necessary to calculate a higher order correction 
to the Hamiltonian. In addition, the conservation properties derived so far did not consider 
collisions. Collisions are important for the radial flux of momentum transport even for very 
low coUisionality because there is a piece of Fi2 that is turbulent in nature and depends on 
the collision operator [35|, |5J]. This piece will give a contribution to Hq in (TT6|) as large as 
any of the other turbulent contributions. The treatment of the collision operator in full / 
simulations remains a challenge. 



IV. DISCUSSION 

It is clear that to reproduce the variety of behaviors observed in intrinsic rotation exper- 
iments with gyrokinetic simulations, it is necessary to retain many effects that are small in 
p=K compared to the lowest order gyrokinetic formulation. These effects that are in principle 
small in p* become the dominant contribution because of a cancellation of the lowest order 
momentum flux due to a symmetry in the system. Many new contributions that have been 
previously neglected now become important. To attain predictive capability, we need to let 
the different contributions compete to determine the rotation direction. The problem is very 
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complex, but it can be simplified using a subsidiary expansion in Bp/ B <^ 1 [35|, l5J]. This 
expansion has allowed us to prove that the intrinsic rotation will depend on the density and 
temperature gradients 5J] , and it has given us a characteristic size of the intrinsic rotation 
generated by the turbulence, Vj ~ [B / Bp)p^Vti- This order of magnitude estimate gives a 
clear prediction for the scaling of the intrinsic rotation generated in the core of tokamaks 
with temperature and plasma current that has been confirmed by comparing data from 
several machines [16|. 

The implementation of the new low flow formulation of gyrokinetics necessary to describe 
intrinsic rotation is challenging, but it is tractable for existing 5f simulations. The same 
cannot be said for gyrokinetic simulations that evolve the full distribution function. These 
simulations need a gyrokinetic equation derived to high order in p^,, and they cannot benefit 
from many of the simplifications that make 5f simulations computationally faster and more 
manageable. The required high order gyrokinetic equation can be obtained in the collision- 
less limit if the gyrokinetic Hamiltonian is obtained to order p^Tg. The addition of collisions 
to these full / models remains unresolved. It is very likely that the full / gyrokinetic for- 
mulation can be simplified with a subsidiary expansion in Bp/B ^ 1. However, even then 
the numerical accuracy to which the gyrokinetic equation needs to be solved is extremely 
high, and scales like p^, as indicated by the requirement in (128|) . Numerical methods that 
ensure this accuracy need to be considered. Wherever the scale length L of the phenomena 
of interest is small (e.g. transport barriers), the accuracy constraint p^ = {pi/LY may not 
seem as daunting, but it remains a condition to be dealt with, in particular regarding the 
equations, which need to be derived to this order. 
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